home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / integration / qcheb.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-12-09  |  5.7 KB  |  230 lines
  1. /* integration/qcheb.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <gsl/gsl_math.h>
  23. #include <gsl/gsl_integration.h>
  24.  
  25.  
  26. /* This function computes the 12-th order and 24-th order Chebyshev
  27.    approximations to f(x) on [a,b] */
  28.  
  29. void
  30. gsl_integration_qcheb (gsl_function * f, double a, double b, double *cheb12, double *cheb24)
  31. {
  32.   size_t i;
  33.   double fval[25], v[12];
  34.  
  35.   /* These are the values of cos(pi*k/24) for k=1..11 needed for the
  36.      Chebyshev expansion of f(x) */
  37.   
  38.   const double x[11] = { 0.9914448613738104,     
  39.              0.9659258262890683,
  40.              0.9238795325112868,     
  41.              0.8660254037844386,
  42.              0.7933533402912352,     
  43.              0.7071067811865475,
  44.              0.6087614290087206,     
  45.              0.5000000000000000,
  46.              0.3826834323650898,     
  47.              0.2588190451025208,
  48.              0.1305261922200516 };
  49.   
  50.   const double center = 0.5 * (b + a);
  51.   const double half_length =  0.5 * (b - a);
  52.   
  53.   fval[0] = 0.5 * GSL_FN_EVAL (f, center + half_length);
  54.   fval[12] = GSL_FN_EVAL (f, center);
  55.   fval[24] = 0.5 * GSL_FN_EVAL (f, center - half_length);
  56.  
  57.   for (i = 1; i < 12; i++)
  58.     {
  59.       const size_t j = 24 - i;
  60.       const double u = half_length * x[i-1];
  61.       fval[i] = GSL_FN_EVAL(f, center + u);
  62.       fval[j] = GSL_FN_EVAL(f, center - u);
  63.     }
  64.  
  65.   for (i = 0; i < 12; i++)
  66.     {
  67.       const size_t j = 24 - i;
  68.       v[i] = fval[i] - fval[j];
  69.       fval[i] = fval[i] + fval[j];
  70.     }
  71.  
  72.   {
  73.     const double alam1 = v[0] - v[8];
  74.     const double alam2 = x[5] * (v[2] - v[6] - v[10]);
  75.  
  76.     cheb12[3] = alam1 + alam2;
  77.     cheb12[9] = alam1 - alam2;
  78.   }
  79.  
  80.   {
  81.     const double alam1 = v[1] - v[7] - v[9];
  82.     const double alam2 = v[3] - v[5] - v[11];
  83.     {
  84.       const double alam = x[2] * alam1 + x[8] * alam2;
  85.  
  86.       cheb24[3] = cheb12[3] + alam;
  87.       cheb24[21] = cheb12[3] - alam;
  88.     }
  89.  
  90.     {
  91.       const double alam = x[8] * alam1 - x[2] * alam2;
  92.       cheb24[9] = cheb12[9] + alam;
  93.       cheb24[15] = cheb12[9] - alam;
  94.     }
  95.   }
  96.  
  97.   {
  98.     const double part1 = x[3] * v[4];
  99.     const double part2 = x[7] * v[8];
  100.     const double part3 = x[5] * v[6];
  101.     
  102.     {
  103.       const double alam1 = v[0] + part1 + part2;
  104.       const double alam2 = x[1] * v[2] + part3 + x[9] * v[10];
  105.       
  106.       cheb12[1] = alam1 + alam2;
  107.       cheb12[11] = alam1 - alam2;
  108.     }
  109.     
  110.     {
  111.       const double alam1 = v[0] - part1 + part2;
  112.       const double alam2 = x[9] * v[2] - part3 + x[1] * v[10];
  113.       cheb12[5] = alam1 + alam2;
  114.       cheb12[7] = alam1 - alam2;
  115.     }
  116.   }
  117.  
  118.   {
  119.     const double alam = (x[0] * v[1] + x[2] * v[3] + x[4] * v[5]
  120.         + x[6] * v[7] + x[8] * v[9] + x[10] * v[11]);
  121.     cheb24[1] = cheb12[1] + alam;
  122.     cheb24[23] = cheb12[1] - alam;
  123.   }
  124.  
  125.   {
  126.     const double alam = (x[10] * v[1] - x[8] * v[3] + x[6] * v[5] 
  127.         - x[4] * v[7] + x[2] * v[9] - x[0] * v[11]);
  128.     cheb24[11] = cheb12[11] + alam;
  129.     cheb24[13] = cheb12[11] - alam;
  130.   }
  131.  
  132.   {
  133.     const double alam = (x[4] * v[1] - x[8] * v[3] - x[0] * v[5] 
  134.         - x[10] * v[7] + x[2] * v[9] + x[6] * v[11]);
  135.     cheb24[5] = cheb12[5] + alam;
  136.     cheb24[19] = cheb12[5] - alam;
  137.   }
  138.  
  139.   {
  140.     const double alam = (x[6] * v[1] - x[2] * v[3] - x[10] * v[5] 
  141.         + x[0] * v[7] - x[8] * v[9] - x[4] * v[11]);
  142.     cheb24[7] = cheb12[7] + alam;
  143.     cheb24[17] = cheb12[7] - alam;
  144.   }
  145.  
  146.   for (i = 0; i < 6; i++)
  147.     {
  148.       const size_t j = 12 - i;
  149.       v[i] = fval[i] - fval[j];
  150.       fval[i] = fval[i] + fval[j];
  151.     }
  152.  
  153.   {
  154.     const double alam1 = v[0] + x[7] * v[4];
  155.     const double alam2 = x[3] * v[2];
  156.  
  157.     cheb12[2] = alam1 + alam2;
  158.     cheb12[10] = alam1 - alam2;
  159.   }
  160.  
  161.   cheb12[6] = v[0] - v[4];
  162.  
  163.   {
  164.     const double alam = x[1] * v[1] + x[5] * v[3] + x[9] * v[5];
  165.     cheb24[2] = cheb12[2] + alam;
  166.     cheb24[22] = cheb12[2] - alam;
  167.   }
  168.  
  169.   {
  170.     const double alam = x[5] * (v[1] - v[3] - v[5]);
  171.     cheb24[6] = cheb12[6] + alam;
  172.     cheb24[18] = cheb12[6] - alam;
  173.   }
  174.  
  175.   {
  176.     const double alam = x[9] * v[1] - x[5] * v[3] + x[1] * v[5];
  177.     cheb24[10] = cheb12[10] + alam;
  178.     cheb24[14] = cheb12[10] - alam;
  179.   }
  180.  
  181.   for (i = 0; i < 3; i++)
  182.     {
  183.       const size_t j = 6 - i;
  184.       v[i] = fval[i] - fval[j];
  185.       fval[i] = fval[i] + fval[j];
  186.     }
  187.  
  188.   cheb12[4] = v[0] + x[7] * v[2];
  189.   cheb12[8] = fval[0] - x[7] * fval[2];
  190.  
  191.   {
  192.     const double alam = x[3] * v[1];
  193.     cheb24[4] = cheb12[4] + alam;
  194.     cheb24[20] = cheb12[4] - alam;
  195.   }
  196.  
  197.   {
  198.     const double alam = x[7] * fval[1] - fval[3];
  199.     cheb24[8] = cheb12[8] + alam;
  200.     cheb24[16] = cheb12[8] - alam;
  201.   }
  202.  
  203.   cheb12[0] = fval[0] + fval[2];
  204.  
  205.   {
  206.     const double alam = fval[1] + fval[3];
  207.     cheb24[0] = cheb12[0] + alam;
  208.     cheb24[24] = cheb12[0] - alam;
  209.   }
  210.  
  211.   cheb12[12] = v[0] - v[2];
  212.   cheb24[12] = cheb12[12];
  213.  
  214.   for (i = 1; i < 12; i++)
  215.     {
  216.       cheb12[i] *= 1.0 / 6.0;
  217.     }
  218.  
  219.   cheb12[0] *= 1.0 / 12.0;
  220.   cheb12[12] *= 1.0 / 12.0;
  221.  
  222.   for (i = 1; i < 24; i++)
  223.     {
  224.       cheb24[i] *= 1.0 / 12.0;
  225.     }
  226.  
  227.   cheb24[0] *= 1.0 / 24.0;
  228.   cheb24[24] *= 1.0 / 24.0;
  229. }
  230.